*	See Gibbs estimation of microstructure models: Teaching notes
	November, 2006
	Joel Hasbrouck
____________________________________________________________________________________________________

	RollGibbsBeta: macro to estimate Roll model with Gibbs sampler

	RollGibbs(
		dsIn, 				input dataset with price variable 'p'
		qOut=qOut, 			output dataset containing simulated q's
		parmOut=parmOut, 	output dataset containing simulated parameters
		nSweeps=100,		number of sweeps
		qDraw=1, 			Draw (simulate) the q's (set qDraw=0 otherwise)
		varuDraw=1, 		Draw varu
		regDraw=1, 			Draw the regression coefficient (c)
		varuStart=0.001, 	Starting value for varu
		cStart=.01, 		Starting value for c
		printLevel=0)

____________________________________________________________________________________________________;


%let Infinity=1e30;
%let eps=1e-30;

%macro RollGibbsBeta(dsIn, qOut=qOut, parmOut=parmOut, nSweeps=100,
	qDraw=1, varuDraw=1, regDraw=1, varuStart=0.001, cStart=.01, betaStart=1, printLevel=0,
	cLower=0, cUpper=&Infinity);
proc iml;
	start main;
	reset printadv=1;
	call streaminit(1234);	*	Initialize the random number generators;
	reset storage=this.imlstor; load;	*	Reload necessary subroutines;
	show modules;

	*	Read in data;
	use &dsIn;
	if &qDraw=1 then do;
		read all var {p pm};
		q = j(nrow(p),1,.);
	end;
	else read all var {p q pm};
	nObs = nrow(p);
	print p q pm;

	*	Create output datasets for simulated qs and parameters;
	create &qOut var {sweep t q};
	create &parmOut var {sweep sdu c beta};

	dp = p[2:nObs] - p[1:(nObs-1)];

	if &qDraw then do; 	*	Initialize qs to price sign changes;
		qInitial = {1} // sign(dp);
		qInitial = qInitial # (q^=0);	*	Only initialize nonzero elements of q;
		q = qInitial;
	end;

	varu = &varuStart;
	c = &cStart;
	beta = &betaStart;

	do sweep=1 to &nSweeps;

		if mod(sweep,500)=0 then print sweep [label='' rowname='Sweep:'];

		dq =  q[2:nObs] - q[1:(nObs-1)];
		dpm = pm[2:nObs] - pm[1:(nObs-1)];

		if &regDraw then do;
			priorMu = {0,1};
			postMu = priorMu;
			priorCov = diag({1,2});
			postCov = priorCov;
			X = colvec(dq) || colvec(dpm);
			rc = BayesRegressionUpdate(priorMu, priorCov, dp, X, varu, postMu, postCov);
			if &printLevel>=2 then print postMu postCov;
			coeffLower={&cLower,-&Infinity};
			coeffUpper=j(2,1,&Infinity);
			coeffDraw = mvnrndT(postMu, postCov, coeffLower, coeffUpper);
			if &printLevel>=2 then print coeffDraw;
			c = coeffDraw[1];
			beta = coeffDraw[2];
		end;

		if &varuDraw then do;
			u = dp - c*dq - beta*dpm;
			priorAlpha = 1.e-12;
			priorBeta = 1.e-12;
			postAlpha = .;
			postBeta = .;
			rc = BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta);
			x = (1/postBeta) * rand('gamma',postAlpha);
			varu = 1/x;
			sdu = sqrt(varu);
			if &printLevel>=2 then print varu;
		end;

		if &qDraw then do;
			qDrawPrintLevel = 0;
			p2 = p - beta*pm;
			call qDraw(p2, q, c, varu, qDrawPrintLevel);
			setout &qOut;
			mOut = j(nObs,1,sweep) || colvec(1:nObs) || colvec(q);
			append from mOut;
		end;

		if &regDraw | &varuDraw then do;
			setout &parmOut;
			append var _all_;
		end;
	end;
	finish main;

run;
quit;
%mend RollGibbsBeta;
